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Abstract 

This  paper  investigates  the  channel  flow  with  a  suction  or  injection  boundary  condition  to  imitate  the  reactant  flow  in  fuel  cells.  A  two- 
dimensional  model  is  proposed  for  a  channel  with  a  trapezoidal  cross-section,  and  systematically  analyzed  using  the  perturbation  method  and  the 
Laplace  transform  technique.  The  model  predicts  the  species  concentration  at  the  electrode  of  the  fuel  cell  under  the  constraint  of  uniform  current 
density  along  the  channel.  For  convenience  of  use  as  a  formula,  the  asymptotic  solutions  to  some  results  are  found  for  two  different  limits  of  xlk. 
©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Channel  flow  belongs  to  the  class  of  most  fundamental 
problems  for  which  the  Navier-Stokes  equation  can  be  solved 
exactly.  Since  every  device  using  fluids  needs  a  fluid  supply 
passage,  channel  flow  has  functioned  as  a  simplistic  model  of 
the  phenomena  that  occur  in  complex  devices  such  as  catalytic 
converter-fuel  cell  assemblies  and  MEMS  devices.  Considering 
its  wide  applicability  and  versatility,  it  might  not  be  possible  to 
enumerate  all  the  applications  of  channel  flow.  As  models  have 
been  developed  to  incorporate  various  aspects  of  physics,  tech¬ 
niques  to  solve  them  with  greater  sophistication  and  speed  have 
also  been  suggested.  However,  in  spite  of  the  long  history  of 
development,  most  of  these  models  do  not  allow  an  analytical 
approach  because  of  the  non-linearity  of  the  governing  equation 
and  the  coupling  between  variables. 

A  fuel  cell  is  an  electrochemical  energy-conversion  device 
that  uses  different  types  of  fuel,  such  as  hydrocarbons,  methanol 
and  hydrogen.  There  are  various  types  of  fuel  cell,  depending 
on  the  electrolyte,  but  most  of  them  have  a  supply  channel  of 
reactants  to  a  porous  electrode  in  common.  From  the  perspec¬ 
tive  of  the  fuel  cell  channel,  the  electrochemical  reaction  at  the 
electrode  can  be  modeled  as  suction  of  reactants  and  injection  of 
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product  species  through  the  bottom  of  the  channel.  The  channel 
flow  with  wall  suction  or  injection  without  the  inclusion  of  a 
chemical  reaction  is  well  known  as  the  Berman  problem.  Many 
researchers  have  tried  to  find  various  types  of  solutions  and  the 
instability  modes  for  the  so-called  symmetric  and  asymmetric 
Berman  problems.  A  two-dimensional  velocity  profile  has  been 
found  for  the  Berman  problem  with  the  aid  of  a  similarity  func¬ 
tion  that  can  be  expressed  as  a  series  expansion  in  terms  of 
the  Reynolds  number  based  on  the  injection  or  suction  veloc¬ 
ity.  Recently,  Cox  and  King  [1]  reviewed  this  approach  in  their 
findings  of  more  complete  asymptotic  solutions  for  the  Berman 
problem.  In  the  field  of  fuel  cells,  where  electrochemical  reac¬ 
tions  complicate  the  channel  flow,  Kulikovsky  [2]  analytically 
found  the  velocity  and  species  concentration  along  the  channel 
direction.  Using  a  one-dimensional  model  and  neglecting  vis¬ 
cous  dissipation,  he  used  a  low-Mach-number  approximation  to 
remove  the  non-linearity  of  the  convection  term. 

In  this  paper,  we  propose  a  two-dimensional  reacting 
flow  model  for  a  fuel  cell  that  includes  variation  of  species 
concentration  in  the  direction  of  the  channel  height.  The 
channel  is  allowed  to  have  a  trapezoidal  cross-section  instead 
of  a  rectangular  shape  so  that  the  model  can  predict  the  effect  of 
the  draft  angle  that  is  common  in  composite  or  metal  separator 
plates.  In  this  paper  it  is  assumed  that  the  velocity  is  not  changed 
by  the  electrochemical  reaction,  which  can  be  justified  when  the 
concentration  of  inert  gas  is  exceedingly  large  compared  to  that 
of  the  species  participating  in  the  reaction  [2] .  This  hypothesis 
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Nomenclature 

A,  A'  constant  representing  the  concentration  gradient 
c  species  concentration 

C  Laplace  transform  of  c 

F  Faraday’s  constant 

h  channel  half-height 

i  current  density 

j  index  used  in  the  Stehfest  algorithm  of  Eq.  (31) 

k  parameter  defined  in  Section  3.1 

L  Laplace  transform 

M  number  of  terms  in  the  Stehfest  algorithm  of  Eq. 

(31) 

mWaii  mass  flux  at  the  bottom  wall 

O  order  of  magnitude 

p  pressure 

Re  Reynolds  number 

s  Laplace-transformed  coordinate  of  a 

Sc  Schmidt  number 

u  stream- wise  velocity 

v  transverse  velocity 

W  channel  width 

Wo  channel  width  at  y  =  0 

Wm,u  constant  defined  in  Section  4.1 

v  non-dimensionalized  stream-wise  coordinate 

y  non-dimensionalized  transverse  coordinate 

Greek  letters 

(3  draft  angle  parameter  defined  in  Section  2 

s  perturbation  parameter 

A  constant  defined  in  Section  2 

0  draft  angle 

Subscripts 
0,(0)  zero  order 

1, (1)  first  order 

2, (2)  second  order 

n  index  used  in  the  Stehfest  algorithm  of  Eq.  (31) 

w  wall  condition 

Superscripts 
i  species  i 

n  index  used  in  the  Stehfest  algorithm  of  Eq.  (31) 


pvW  (y)  A*  +  pvW  (y)  Ay+ 

d/dy  [pv\N  (y)  Ax]  Ay  d/dx  [pu  w  (y)  Ay]  Ax 


Fig.  2.  Infinitesimal  control  volume  for  the  channel  with  draft  angle  6. 


needs  to  be  checked  for  many  examples  of  real  operating 
conditions,  but  is  adopted  for  simplicity  of  the  solution  in  this 
paper. 


2.  Mathematical  formulation 


As  a  method  to  investigate  the  reacting  flow  in  a  fuel  cell,  we 
use  channel  flow  with  a  suction  or  injection  boundary  condition 
for  a  certain  species.  Fig.  1  shows  a  schematic  diagram  of  the 
model  that  is  investigated  in  this  paper.  The  flow  is  driven  by  a 
pressure  gradient  formed  in  the  a  direction,  while  a  concentration 
gradient  in  the  y  direction  imparts  additional  movement  to  the 
species  within  the  channel.  In  the  z  direction,  there  are  variations 
in  the  velocity  and  concentration  originating  from  the  no- slip  and 
zero-flux  boundary  condition  at  the  wall.  However,  its  effect  is 
secondary  compared  with  variation  in  the  a  and  y  directions. 
The  fuel  cell  channel  has  a  small  cross-section  compared  with 
its  length,  which  often  amounts  to  a  few  100  times  its  cross- 
section  dimension.  The  channel  width  W(y)  varies  with  y  and 
can  be  written  as  function  of  the  draft  angle  0 ,  for  which  0  =  0 
corresponds  to  the  rectangular  channel  case. 

In  the  fuel  cell,  the  suction  velocity  of  the  reactant  species 
is  a  function  of  the  current  density  divided  by  the  product  of 
the  Faraday  constant  and  the  concentration  of  that  species  [3]. 
Considering  that  the  average  current  density  of  a  fuel  cell  under 
normal  operating  conditions  is  approximately  1  A  cm-2,  the  suc¬ 
tion  velocity  is  found  to  be  a  tiny  fraction  of  the  mainstream 
velocity  U,  which  is  approximately  5  m  s-1  under  the  same  con¬ 
ditions. 


F  w 

77 


o(io-3), 


3  d 
3 y  ^  3a 


This  observation  allows  us  to  assume  that  diffusion  in  the  y 
direction  will  dominate  that  in  the  other  directions,  which  is 
consistent  with  the  boundary  layer  approximation  [4].  Fig.  2 
shows  the  control  volume  for  the  trapezoidal  channel  in  Fig.  1, 


Streamwise  Concentration 
Velocity  Profile  Profile 


Fig.  1.  Schematic  diagram  of  the  channel  shape  and  coordinate  system  under  investigation. 
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together  with  mass  fluxes  into  and  out  of  each  surface.  It  has 
finite  width  W(y ),  and  infinitesimal  thickness  Ax  and  height  Ay. 
Taking  the  limit  as  Ax  — ^  0,  Ay  — >  0,  the  mass  balance  within 
the  control  volume  yields: 


3(Wu)  +  d(Wv)  _  Q 


3x 


dy 


The  momentum  and  the  species  balance  can  be  treated  in  the 
same  manner.  As  a  first  approximation  to  the  Maxwell-Stefan 
equation,  the  diffusion  is  considered  to  follow  Fick’s  law  for 
binary  fluids,  irrespective  of  the  number  of  components.  The 
length,  velocity,  and  pressure  are  non-dimensionalized  by  the 
channel  half-height  h ,  mean  velocity  U,  and  dynamic  pres¬ 
sure  pU2.  Neglecting  axial  diffusion  and  introducing  non- 
dimensionalized  variables,  the  governing  equation  for  the  trape¬ 
zoidal  channel  flow  with  reaction  yields: 


-  du  -  dv  dW 

W —  +  W —  +  V -  =  0 

3x  dy  dy 


_  du 

-  du 

-  dp 

1  d 

Wu  — 

+  Wn—  = 

W  y  + 

— 

9x 

dy 

3x 

Re  dy 

_  dcl 

-  dcl 

1  d 

-  dcl 

Wu 

+  Wp 

W 

dx 

dy 

Re  Sc  dy 

9y  _ 

-  du 


The  normalized  channel  width  W(y)  and  the  draft  angle  param¬ 
eter  p  are  defined  as  follows: 

2  tan  0 

W(y)=\-py  (-1  <  y  <  1),  P  =  r 

Wo/  h 

where  6  is  the  draft  angle,  and  Wo  =  W (0)  (4) 


Since  we  assume  that  the  velocity  is  not  changed  by  the  reaction, 
the  momentum  equation  for  the  fully  developed  flow  can  be 
simplified  as  follows: 


du 

dx 


v  =  0  •  •  •  Re 


dp 

dx 


=  const  =  X 

(5) 


From  the  no-slip  boundary  condition: 


u(- 1)  =  u(  1)  =  0 


Expanding  u  in  terms  of  the  draft  angle  parameter  ft  as  a  pertur¬ 
bation  variable: 


u  —  u(0)  +  Pu(l)  +  /32W(2)  +  •  •  • , 

where  the  zero  order  in  is: 
d  U(0) 

^  2  =  w(0)(-l)  =  «(0)(1)  =  0 

X  2 

M(0)  =  2  (y  - 


the  first  order  in  fi  is: 


U(i)(—1)  =  U(  i)(l)  =  0 

the  second  order  in  ft  is: 

d2«(2)  =  _d_  /  dw(i)\ 
dy2  dy  V  dy  ) 


K(2)(  — 1)  =  M(2)(l)  =  0 


Ay  =  Ay, 


M(i)  =  ^(y3  -  y)  (9) 

X  o 

3  (2y3  -  y), 

“(2)  =  ^3;y5  “  5y3  +  2y ^ 

(10) 


and  the  third  order  in  p  is: 


d2^(3) 


y 


dM(2)\ 


=  T(75y4  -  45 y 2  +  2), 


dy2  dy 
W(3)(— 1)  =  «(3)(1)  =  0 

•  ■.«(3)  =  A(10/-15y4  +  4y2+1) 


(ID 


Collecting  and  adding  terms  up  to  the  third  order  gives: 

(3y5  -  5y3  +  2y) 


2  i  \  i  q/Xj  s  3  \  I  £>2  ^  /o  5  r  3 


U  K  ^(y2-  l)  +  AyJ  -y)  +  ,. 

2  6  90 


+p- 


A 


360 


(10y6  -  15y4  +  4y2  +  1) 


(12) 


The  constant  A  is  calculated  from  the  definition  of  the  average 
velocity: 


X  ^  —3  (13) 


Fig.  3  shows  the  change  in  velocity  profile  with  draft  angle  var¬ 
ied  up  to  40°  for  the  case  Wo/h  =  2.  The  location  of  the  peak 
velocity  moves  upward  as  the  draft  angle  increases  because  of 
the  reduction  in  the  effective  area  in  the  upper  region  of  channel. 


1) 


(8) 


Fig.  3.  Effect  of  draft  angle  6  on  normalized  velocity  profile  within  the  channel. 


1054 


Y.  Sung  /  Journal  of  Power  Sources  159  (2006)  1051-1060 


This  is  the  result  for  Wo/h  =  2,  in  other  words  for  a  channel  with 
square  cross-section  when  it  has  no  draft  angle.  If  the  channel 
has  greater  width  than  its  height,  the  effect  of  the  draft  angle  will 
be  lessened,  as  can  be  deduced  from  the  definition  of  /3. 

Using  Eqs.  (12)  and  (13),  the  species  equation  now  can  be 
simplified  as  follows: 


_  3  cl 

Wu  — 
dx 


1  3 

Re  Sc  3 y 


-  3 cl 
W  — 


where  u  &  U(p)  +  /*«(i)  +  P2^(2)  +  /33U( 3)  H - 


(14) 


As  for  the  boundary  condition,  the  inlet  concentration  is  assumed 
to  be  constant.  To  simplify  the  calculation,  we  redefine  the  con¬ 
centration  cl  as  deviation  from  the  inlet  condition.  In  spite  of  this 
change,  Eq.  (14)  can  still  be  used  in  its  original  form,  since  it  con¬ 
tains  only  derivative  terms.  At  the  reacting  wall  of  the  bottom, 
the  Neumann-type  boundary  condition  can  be  imposed  under  the 
assumption  that  mass  transfer  there  occurs  only  by  diffusion.  If 
the  channel  has  a  meandering  shape  or  a  dead-end  configura¬ 
tion,  the  pressure  rise  in  the  channel  will  be  high  enough  to 
allow  the  direct  penetration  of  species  into  the  electrode.  In  that 
case,  we  have  to  include  convective  transport  into  the  wall,  which 
requires  simultaneous  solution  of  the  y-momentum  equation.  To 
simplify  the  solution,  we  assume  a  straight  channel  that  has  neg¬ 
ligible  pressure  drop  along  the  channel.  Since  the  total  flux  of 
suction  or  injection  is  determined  by  external  loading  and  is  not 
dependent  on  the  change  in  the  bottom  area  of  channel  caused  by 
draft,  the  following  condition  is  imposed  on  the  bottom  surface: 


(d  cl 

Wwall  ~  Wy=-l  — 

V  fy 


A 

1+/3 


where  A  ~  mwa n 


(15) 


At  the  upper,  non-reacting  wall,  there  will  be  no  mass  flux  of 
any  kind,  which  means  that  the  concentration  gradient  vanishes 
there.  To  sum  up: 


3  cl 

—  =  0 
fy  y=l 


(16) 


If  A  is  positive,  suction  of  species  i  occurs  at  the  wall,  whereas 
if  it  is  negative,  injection  occurs.  The  magnitude  of  A  is  related 
to  the  amount  of  suction  or  injection  of  species  i  and  conse¬ 
quently  to  the  local  current  density  by  Faraday’s  law.  Under  most 
operating  conditions  for  fuel  cells,  except  at  very  low  current 
density,  A  varies  along  x  and  should  be  simultaneously  deter¬ 
mined  with  the  other  variables.  Although  the  current  density  is 
sometimes  assumed  to  have  a  certain  profile,  such  as  an  expo¬ 
nentially  decreasing  function  [2],  or  to  have  constant  value  [5], 
it  is  obvious  that  assumption  of  the  condition  of  constant  A  is 
not  well  supported  by  experiments  in  most  cases.  However,  we 
adopt  it  in  this  paper  to  concentrate  on  the  flow  and  mass  transfer 
while  considering  the  electrochemical  reaction  only  through  the 
suction  or  injection  boundary  condition. 

The  concentration  of  species  /,  more  exactly  the  concentration 
change  for  species  i  from  the  inlet  condition,  can  be  expanded 


in  terms  of  ft  as  in  the  case  of  velocity: 


°l  ~  c(0)  +  Pc\\)  +  ^c\l)  +  P3 c(3)  + 


(17) 


Expanding  Eqs.  (14)  and  (16)  similarly  and  collecting  like 
terms  of  ft,  the  following  equations  and  boundary  conditions 
are  obtained  after  substitution  into  Eqs.  (8)-(l  1). 

Zero-order  terms  in  ft: 


Kl-y2)- ^ 

dx 


2J 


3  c 


(0) 


3  r 


where  k  =  -Re  Sc 
2 


(18) 


c 


(0) 


x=0 


=  0, 


9c 


(0) 


3  y 


A 


y=- 1 


1  + 


3c 


(0) 


9  y 


=  0 


y=i 


(19) 


First-order  terms  in  ft: 


c 


(i) 


A'=0 


x  9c(l) 

;  dx 

d2cl 

0  c(1) 

dy2 

_  9c(0) 
dy 

1  „a  c(0) 
3 37  dy2 

(20) 

=  0. 

9cU) 

dy 

=  0, 

y=- 1 

dcW 

dy 

=  0 

y=l 

(21) 

Second-order  terms  in  /3: 


kd-y2)^ 

dx 


2  J 


d  c 


(2) 


3c 


(i) 


y  aMi) 


2y  3c 


(0) 


dy  3  dyx 


2y  9  9  C/q'j 

+  — (2  +  5y-3y2) - ^ 

90  *  *  >  dy2 


3  dy 


(22) 


c 


(2) 


x=0 


=  0, 


3c 


(2) 


dy 


=  0, 


y=-i 


Third-order  terms  in  ft: 


k(l  - 


dr1 

Pit 

dx 


2J 


3  c 


(3) 


3c 


(2) 


dy2 

y  ^c(2) 


3  ci 


(2) 


dy 


=  0  (23) 


y=i 


3  y  3  3 y 


2y  9cqx  y 

—  — ^  +  ^-(2  +  5y 
3  3y  45  ^ 


d2cl 

3r)^f 

dy1 


y  „  -  2^9c(o> 

—  tt(2  +  —  3y  )— — 

45  dy 

+^(3~y2- 2  °y3  -  6V) 


d2cl 

4^d  C(0) 
dy2 


(24) 


c 


(3) 


x=0 


=  0, 


3c 


(3) 


dy 


=  0, 


9c 


y=-i 


(3) 


dy 


=  0  (25) 


y=i 


3.  Numerical  results 


3.1.  The  perturbation  method 


In  this  section,  we  use  the  perturbation  method  and  the 
Laplace  transform  technique  to  systematically  investigate  Eqs. 
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(1 8)— (25).  The  perturbation  method  in  this  section  is  different 
from  that  used  in  Section  2,  in  that  the  perturbation  parameter 
s  has  no  physical  meaning,  but  is  introduced  for  mathematical 
convenience.  Eqs.  (18),  (20),  (22)  and  (24)  commonly  contain 
the  term  y 2  on  the  left-hand  side,  the  omission  of  which  greatly 
simplifies  the  solution.  Hence,  we  multiply  the  term  y2  by  the 
perturbation  parameter  e  and  expand  the  left-hand  side  in  terms 
of  £,  while  the  right-hand  side  is  calculated  from  the  solutions 
for  the  lower-order  equations.  After  solving  the  decomposed 
equations,  they  are  recombined  with  the  substitution  of  1  for 
the  perturbation  parameter  e.  This  approach  can  be  regarded  as 
being  valid  as  far  as  the  recombined  terms  constitute  a  converged 
series  [6,7]. 

3.1.1.  Zero-order  terms  in  /3 


kil-sy2)^ 

dx 


2  J 


dZC 


(0) 


dy * 


where  C(0)  =  co  +  sc\  +  £  C2  + 

Zero-order  terms  in  £\ 


(26) 


3co  _  d2co 
"to  ~  ~dy2' 


CQ  lx=0  0, 


dcp 

dy 


A 


y=- 1 


1  + 


dco 

dy 


=  0 


(27) 


y=i 


For  terms  of  order  j  larger  than  0: 

dcj- 1 


dcj  d2Cj  , 
k — - — — r=kyx 


dx  d y 


dx 


CJ 


x=0 


dc  j 

o,  — 

dy 


=  0, 


y=—i 


dcj_ 

dy 


=  0  where  j  =  1,  2,  3, 


(28) 


y=i 


3.1.2.  First-order  terms  in  j3 


?  dcU) 

H  l  -  ey2)-^ 
dx 


2  J 


d  c 


(1) 


del 


(0) 


y  aMo) 


dy< 


dy  3  dyx 


where  c1^  =  cq  +  £C\  +  £2C2  + 
Zero-order  terms  in  £\ 


(29) 


3co  a2co 

k - 7T 

dx  dy2 


dc 


(0) 


dco 

dy 


=  0, 


_ y  8Mo) 

dy  3  dy2 
dco 


^0  I  r=0  6. 


=  0 


y=i 


y=—l 

For  terms  of  order  j  larger  than  0: 

dCj- 1 


(30) 


dcj  d2Cj 

k — - 3r 

dx  dy2 


=  ky 


dx 


ci 


_n  =  0, 


x=0 


dc 


J 


dy 


=  0, 


y=- 1 


dcj 

dy 


=  0  where  j  =  1,  2,  3, 


(31) 


y=i 


3.1.3.  Second-order  terms  in  /3 


k{\  -cy2)^ 
dx 


2  J 


dzc 


(2) 


dc 


(1) 


y  aMi) 


2y 


(0) 


dy< 


dy  3  dy* 


3  9y 


2y  9  d2  ci o) 

H-  — (2  +  5y  —  3y2) - 22 

90  y  7  dy2 

where  c[2)  =  co  +  sc\  +  e  C2  + 


Zero-order  terms  in  £\ 


(32) 


dco  d2co  dc(i)  y  92C(1)  2y  dc[ 


(0) 


dx 


dy: 


2  y 

+  —  (2  +  5  y 
90  ^ 


dy 


-3  y2) 


3  8y2 

92g(0) 

dy2 


3  dy 


CO  lx=0  0, 


dco 

dy 


=  0, 


y=-i 


dco 

dy 


=  0 


(33) 


y=i 


For  terms  of  order  j  larger  than  0: 


dc  j 
k-1- 
dx 

dcj 


d2C 


dy 


dy2 
=  0, 


=  ky 


dc 


7-1 


dx 

dc 


CJ 


-n=0, 


x=0 


J 


=-l 


y= 


dy 


=  0  where  j  =  1,  2,  3, 


y=i 


(34) 


3.1.4.  Third-order  terms  in  /3 


k{\  -£y2)^ 

av 


2J 


azc 


(3) 


dc 


(2) 


y  a2c[2) 


2y 


(i) 


dy< 


y(2  +  5y-  3y2)  d2c{ 


+ 


(1) 


dy  3  dy2  3  dy 
y( 2  +  35y  -  3y2)  dcl(0) 


45  9j2 

3  -  y2  -  20 y3  -  6 y4  d2cl(0} 
540  dy2 

where  c1^  =  co  +  £C\  +  £  C2  + 


45 


dy 


+ 


(35) 


Zero-order  terms  in  £\ 


k 


dco  92<3Q 


dc\ 


(2) 


y  d2Cl{2) 


2 y 


(i) 


dx 


dy: 


dy  3  dyx 


3  ay 


y(2  +  5y  -  3y2)  92+  _  y(2  +  35y  -  3  y2)  dc 
45  dy2 

3  —  y2  —  20y3  -  6 y4  d2cl(0) 


+ 


+ 

aco 

dy 


(0) 


45 


dy 


540 


dy: 


CO  lx=0  0, 


=  0, 


y=-l 


dco 

dy 


=  0 


(36) 


y=i 
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For  terms  of  order  j  larger  than  0: 


k 


dcj_ 

dx 


d2Cj 


x=0 


dcj_ 

dy 


y=  i 


0  where  j  =  1,2,3,*** 


(37) 


3.2.  The  Laplace  transform 


The  partial  differential  equations  that  were  derived  in  the 
previous  section  can  be  transformed  to  ordinary  differential 
equations  (ODEs)  using  the  Laplace  transform,  the  definition 
of  which  can  be  expressed  as  follows: 


Cj(s,  y)  =  L[cj(x,  y) ]  =  /  e~sxCj(x,  y) dx, 

JO 

where  j  =  0,1,2,*** 


(38) 


3.2.1.  Zero-order  terms  in  f 

C(o)  =  L\c^]  =  Co  +  sC  i  +  s2C2  + 

Zero-order  terms  in  s: 

i2, 


ksCo(s,  y )  = 
dC0 


d  y 


d  CoCf  y) 

dy2 

A  1 


y=-i 


l  +  Bs’ 


dCo 
d  y 


=  0 


y=i 


For  terms  of  order  j  larger  than  0 

d2Cj(s,  y) 


ksC )(s,  y)  - 


dy 


=  ksy  Cj-\(s,  y), 


dCj 

dy 


=  0, 


dC 


j 


y=- 1 


dy 


=  0 


y=i 


(39) 


(40) 


(41) 


3.2.2.  First-order  terms  in  f 

C(i)  =  L[c =  Co  +  eCx  +  £2C2  + 
Zero-order  terms  in  £\ 

d2C0  dC(0)  v  d2C(0) 


ksC  o  — 
dC0 


d  y 


dy2 
=  0, 


(Q) _ 

dy  3  dy2 

dC0 


y=-i 


dy 


=  0 


y=i 


For  terms  of  order  j  larger  than  0: 
d2C j(s,  y) 


ksC j(s,  y)  - 
dC , 


dy 


=  ksy  Cj-i(s,  y). 


dy 


=  0, 


dC 


j 


y=- 1 


dy 


=  0. 


y=i 


(42) 


(43) 


(44) 


3.2.3.  Second-order  terms  in  f 


C( 2)  =  L\cl( 2)]  —  Co  +  eCl  £^C 2  + 


Zero-order  terms  in  e: 


ksCo  — 


d2C0  dC(i)  y  d2C(i)  2>-dC(o) 


dy* 


2  y 

+— (2  +  5y 
90 


3y2) 


dy  3  dy2 
d2C( 


3  dy 


(0) 


dC0 

dy 


=  0, 


dy2 

dC0 


y=-i 


dy 


=  0 


y=i 


For  terms  of  order  j  larger  than  0: 


ksC j(s ,  y) 
dC 


d  CjG  y)  ,  ^ 

- 2 —  =  fay  C/_i(s,  y), 


dy' 


./ 


dy 


=  0, 


y=-i 


dy 


=  0 


y=i 


(45) 


(46) 


(47) 


3.2.4.  Third-order  terms  in  /3 

C(  3)  =  L[c^]  =  Cq  +  eC]  +  e2C2  +  •  •  • 


Zero-order  terms  in  e: 


^Cq  — 


d2C0  dC(2)  y  d2C(2)  2ydC(i) 


(48) 


+ 


+ 


dy2  dy  3  dy2  3  dy 

y(2  +  5y  -  3y2)  d2C(i)  y(2  +  35y  -  3y2)  dC(0) 


45  dy2 

3  -  y2  -  20y3  -  6y4  d2C(0) 


45 


dy 


dCo 

dy 


540 

dC0 

=  °-  -f- 

v— —  I  dy 


dy' 


=  0 


(49) 


y=i 


For  terms  of  order  j  larger  than  0: 


faC j(s,  y)  - 


d 2Cj(s,  y) 


dy* 


=  ksy2Cj-i(s,  y), 


dq 

dy 


=  0, 


>■=—  1 


dC^ 

dy 


=  0 


(50) 


y=l 


3.3.  Solution  of  ODEs  and  the  inverse  Laplace  transform 

The  ODEs  derived  in  the  previous  section  can  be  solved 
analytically  by  successive  substitution.  However,  the  inverse 
Laplace  transforms  for  these  are  nearly  impossible  to  find 
directly  using  the  conversion  formula.  The  Stehfest  algorithm 
[8,9]  is  an  effective  method  to  find  the  inverse  Laplace  trans- 
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Fig.  4.  Change  in  C[0)(x,  1)  as  the  order  of  approximation  in  s  is  changed  (A  =  10, 
&  =  200,  f>-  0).  Exact  solutions  to  ODEs  for  Cj  were  used  in  the  inverse  Laplace 
transform  with  M-  6  for  the  Stehfest  algorithm. 


forms  approximately: 


x 


In  2  In  2 

q(x,  y )  ^  - 2_^WjCi(sj,  y )  whereby  =  j 

x  j= i 

min(j,M/2) 

f  =  (_l)M/2+;  £ 


w 


(2  n)\nM/2 


n=(l+j)/2 


(M/2  -  n)\n\(n  -  1  )!(j  -  n)!(2w  -  j)! 


The  number  of  terms  M  in  the  Stehfest  algorithm  is  usually  set 
higher  than  or  equal  to  six  [10].  In  this  section,  we  provide  some 
results  using  the  internal  functions  of  Mathematica  to  treat  the 
ODEs  and  the  additional  package  developed  by  Cheng  et  al.  [1 1] 
in  implementing  the  Stehfest  algorithm.  The  number  of  terms 
to  be  included  in  the  expansion  for  the  desired  accuracy  is  also 
checked  before  we  try  to  find  the  analytical  solution. 

Let  us  consider  a  polymer  electrolyte  membrane  (PEM)  fuel 
cell  under  normal  operating  conditions  of  1  atm  and  80  °C. 
The  fuel  cell  is  assumed  to  operate  under  uniform  current  den¬ 
sity  mode  at  approximately  1  A  cm-2.  The  stoichiometric  ratio, 
which  represents  how  many  times  reactants  are  supplied  through 
the  inlet  compared  to  the  quantity  that  actually  participates  in 
the  reaction,  is  ordinarily  set  higher  than  or  equal  to  2.0.  The 
channel  cross-section  within  the  separator  plate  is  considered  to 
have  a  square  shape  with  height  and  width  equal  to  2  mm,  and  its 
length  to  be  much  greater,  as  for  most  serpentine  channel  shapes. 
Let  us  assume  that  the  ratio  of  the  channel  length  to  its  height  is 
250,  in  other  words  the  channel  length  non-dimensionalized  by 
the  channel  half-height  h  is  500.  For  this  channel  shape  and  these 
operating  conditions,  the  Reynolds  number  is  approximately 
200.  The  Schmidt  number  for  different  gas  pairs  is  0.56-0.74 
for  N2,  O2,  and  H2O  compounds  [12].  Hence,  k  in  Eq.  (18) 
is  approximately  200-260,  while  x  varies  from  0  to  500.  Fig.  4 
presents  the  concentration  profile  along  the  channel  in  the  case  of 
no  draft  as  the  order  of  approximation  in  s  is  changed.  In  the  case 
of  no  draft,  the  zero-order  term  in  /3,  C|0)  represents  cl.  Fig.  4 

shows  that  C ^  converges  with  the  inclusion  of  higher-order 


terms  in  s,  and  for  the  first-order  approximation  it  approaches 
the  converged  solution  with  8%  error  on  average.  For  the  second- 
order  approximation  in  s,  a  2%  error  applies  to  CJ0)  on  average. 
From  this  result,  we  choose  to  cut  off  the  terms  higher  than  the 
second  order  in  s  in  this  paper.  Fig.  5  presents  the  concentra¬ 
tion  profile  for  the  channel  with  non-zero  draft  angle  with  the 
order  of  approximation  in  /3  changed.  The  result  presented  is 
for  conditions  of  A  =  10,  k-  200,  and  0=10°.  It  shows  that  the 
first-order  approximation  in  /3  is  sufficient  to  provide  an  accu¬ 
rate  profile  of  cl  in  the  case  of  non-zero  draft.  Fig.  6  shows 
the  error  of  the  first-order  approximation  compared  to  the  third- 
order  approximation  at  various  locations  of  the  channel.  Within 
3%  across  the  channel  length,  the  first-order  approximation  in  /3 
provides  a  good  representation  of  higher-order  approximations 
to  the  converged  solution.  Fig.  7  plots  the  concentration  change 
for  species  i  from  the  inlet  condition  in  the  case  of  no  draft.  The 
monitoring  point  is  at  the  reaction  site  located  at  y  =  —  1  along 
the  channel  width  and  k  is  varied  from  100  to  500.  Excluding  the 
rapid  decay  near  the  inlet,  the  concentration  decreases  with  con¬ 
stant  slope  along  the  channel  in  the  case  of  suction.  The  slope 
of  the  decrease  becomes  less  steep  with  higher  k,  the  product 


(51) 


of  the  Reynolds  number  and  the  Schmidt  number.  For  injection 
of  species  i,  its  concentration  would  increase  along  the  channel 
and  the  slope  of  the  profile  would  decrease  as  k  increases.  It 
can  be  deduced  from  the  graph  that  both  increased  inertia  and 
enhanced  mass  diffusion  can  be  used  for  greater  compensation 
of  the  concentration  change  at  the  reaction  site.  Fig.  8  shows 
the  change  in  concentration  profile  within  the  channel  as  k  is 
increased.  The  graph  shows  obvious  variation  in  its  magnitude 


Nondimensionalized  Streamwise  Coordinate,  x 

Fig.  5.  Change  in  cl(x- 1)  as  the  order  of  approximation  in  f  is  changed  (A  =  10, 
k  =  200, 0=10°).  Exact  solutions  to  ODEs  for  Cj  were  used  in  the  inverse  Laplace 
transform  with  M-  6  for  Stehfest  algorithm. 
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Fig.  6.  First-order  approximation  error  in  f  at  various  locations  of  the  reaction 
site. 


Fig.  7.  Concentration  change  for  species  i  from  the  inlet  condition,  cl(x- 1),  as 
k  is  increased  (A  =  10,  f  =  0). 


Fig.  8.  Concentration  change  for  species  i  from  the  inlet  condition,  cz(100,_y),  as 
k  is  increased  (A  =  10,  /3  =  0). 


Fig.  9.  Effect  of  draft  angle  6  on  the  distribution  of  concentration  cl(x- 1)  in  the 
stream- wise  direction  (Wo/h  =  2,  A  =  10,  &  =  200). 


in  the  direction  of  the  channel  height,  which  justifies  the  use  of 
a  two-dimensional  model  for  concentration.  For  k  above  300, 
the  upper  wall  is  hardly  affected  by  reaction  at  the  bottom.  This 
observation  suggests  that  the  measurement  of  species  concen¬ 
tration  through  the  slot  of  the  upper  wall  may  not  be  sufficient 
to  deduce  the  condition  at  the  bottom.  Figs.  9  and  10  show  the 
draft  angle  effect  on  the  concentration  cl  in  the  directions  of  the 
channel  length  and  height.  Fig.  1 1  presents  the  contour  plot  for 
cl  within  the  channel  for  two  different  draft  angles.  The  bigger 
draft  angle  produces  smaller  perturbation  of  the  concentration 
at  the  bottom,  which  means  the  higher  concentration  of  reac¬ 
tant  and  the  faster  removal  of  product  from  the  reaction  site. 
In  the  fuel  cell,  the  enhanced  activity  of  reactant  oxygen  surely 
contributes  to  the  increase  in  performance.  However,  the  fast 
removal  of  product  water  may  cause  membrane  drying,  which 
would  reduce  proton  conductivity  and  lower  the  performance 
when  operated  under  low  humidity  conditions.  In  the  case  of 
high  humidity  conditions,  the  removal  of  product  water  is  favor- 


Fig.  10.  Effect  of  draft  angle  6  on  the  distribution  of  concentration  deviation 
cz(100,_y)  in  the  direction  of  the  channel  height  (Wo/h  =  2,  A  =  10,  k  =  200). 
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Fig.  11.  Concentration  change  for  species  i  from  the  inlet  condition,  cl(x,y), 
within  the  channel  for  different  draft  angles:  (a)  6-  10°  (solid  line);  and  (b) 
6  =  40°  (dashed  line). 


able,  but  if  condensation  starts  to  occur,  the  reduced  momentum 
at  the  bottom  due  to  the  draft  may  not  be  sufficient  to  withdraw 
water  droplets  at  the  interface  with  the  gas  diffusion  layer.  More 
details  on  the  water  problem  for  fuel  cells  can  be  found  in  the 
literature  [13-17]. 


4.  Approximate  analytical  solution 

In  this  section,  we  try  to  find  the  approximate  analytical 
expression  for  the  numerical  results  presented  in  Section  3.  First, 
we  use  two  different  limits  to  approximate  the  ODE  solutions 
and  then  find  their  inverse  Laplace  transforms.  From  the  results 
presented  in  Figs.  5  and  6,  we  determine  to  expand  cl  only  up  to 
the  first  order  in  /3.  The  analytical  solutions  for  ODEs  in  this  sec¬ 
tion  were  found  using  Mathematica  implementing  the  symbolic 
operation. 


4.1.  Asymptotic  solutions  for  ODEs 


4.1.1.  Small  ks  limit 

For  zero-order  terms  in  f ,  C(0)  that  satisfies  Eqs.  (39)-(41) 
has  the  following  approximate  expression  at  y  =  —  1  for  small 
values  of  ks: 


C 


(0) 


Ak 

T+~ 


13 


+ 


46 


18  k2s2  63 ks 


(52) 


Table  1 

WM,n  forM>  16 


77 

Wm,u 

1 

1 

3/2 

0.939437 

2 

0.693147 

5/2 

0.434112 

3 

0.240227 

7/2 

0.120362 

4 

0.0555044 

9/2 

0.0238368 

5 

0.00961786 

11/2 

0.00367086 

6 

0.00133253 

values  of  ks: 

W-I)*4EpSj!.  <*3  =  1.875. 

^  77=3 

~  —1.5,  as  =  1.75,  o' 6  =  —1.375,  =  0.65625 

(54) 

For  first-order  terms  in  ft,  C( i)  that  satisfies  Eqs.  (42)-(44)  has 
the  following  approximate  expression  at  y  =  —  1  for  large  values 
of  ks: 

A  7  12 

Ak  v — >  Otn 

Qi)(s,  -  )  « - >  "  ,  a3  =  0.277344, 

()  1  +  B^kn/2sn/2 

77=3 

ot\  =  1.05859,  as  =  —1.42839,  a 6  =  2.44596, 

a2  =  —4.14128,  a 8  =  7.82894,  ag  =  —13.6265, 

am  =  19.1652,  an  =  -18.6228,  a12  =  9.28491  (55) 


4.2.  Approximate  inverse  Laplace  transform 

Now  we  use  the  Stehfest  algorithm  of  Eq.  (51)  to  find  the 
inverse  Laplace  transforms  for  Eqs.  (52)  and  (53),  and  Eqs.  (54) 
and  (55).  For  convenience,  we  introduce  a  new  constant,  Wm,h- 

m 

EUJ  j 

f  (56> 

7=1  J 

For  M  >  16,  Wm,u  becomes  nearly  independent  of  M,  the  numer¬ 
ical  value  of  which  is  tabulated  in  Table  1  for  different  values  of 
n. 


For  first-order  terms  in  /3,  C(i)  that  satisfies  Eqs.  (42)-(44)  has 
the  following  approximate  expression  at  y  =  —  1  for  small  values 
of  ks: 


4.2.1.  Small  ks  limit  ( large  x/k  limit ) 

This  limit  corresponds  to  the  case  of  large  x/k. 

The  zero-order  terms  in  ft  from  Eqs.  (51)  and  (52)  are: 


Ak  \  13 

1  +/3  |_18 k2s2 


572 

+  2835 ks 


(53) 


4.1.2.  Large  ks  limit 

For  zero-order  terms  in  /3,  C(0)  that  satisfies  Eqs.  (39)-(41) 
has  the  following  approximate  expression  at  y  =  —  1  for  large 


c{0)(x ,  -1)  ^  - 


A 


1  + 


13Wm,2  /*\  46Wm,i 
18 In 2  \k)  +  63 


(57) 


and  the  first-order  terms  in  ft  from  Eqs.  (51)  and  (53)  are: 


c[i) (x,  -1)  «  - 


A 


1  + 


13Wm,2  fx 
18 In 2  \k 


+ 


572  WM,i 
2835 


(58) 
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To  sum  up,  for  large  x/k : 

cl(x,  —  1)  ~  0)0*:,  —1)  of  Eq.  (57)  +  /3cl^(x,  —1)  ofEq.  (58) 

(59) 


c\x,  — 1)  ^  —A 


x 


0.722  -  +  0.730 / (0) 


where  f(/3)  = 


1  +  0.2163/3 


1  + 


(60) 


4.2.2.  Large  ks  limit  (small  x/k  limit ) 

The  zero-order  terms  in  /3  from  Eqs.  (51)  and  (54)  are: 


C(0)(*,  -1)  ~  “ 


A 


7 


£ 


WM,n/20ln  (X\n/ 2-1 


1  +P(f'3  L(ln2)"/2“1  U 


(Z3  =  1.875,  ot\  ~  —1.5,  o' 5  =  1.75,  a 6  =  —1.375, 

0^7  =  0.65625 

and  the  first-order  terms  in  /3  from  Eqs.  (51)  and  (55)  are: 

12 


(61) 


A  ^-^WM,n/2Wn  fX\n/~  1 
ln2)w/2_1  U 


(*3  =  0.277344,  04  =  1.05859,  a 5  =  —1.42839, 

04  =  2.44596,  <27  =  —4.14128,  a%  =  7.82894, 

<29  =  —13.6265,  aqo  =  19.1652,  a\  \  =  —18.6228, 


o'  12  =  9.28491 


(62) 


To  sum  up,  for  small  x/k: 

c\x ,  —  1)  ~  c^(x,  —1)  ofEq.  (61)  +  /3cl^(x,  —1)  ofEq.  (62) 

(63) 


Fig.  12.  Comparison  of  the  asymptotic  approximation  derived  in  Section  4  with 
the  numerical  result  (A  =  10,  k  =  200,  0=10°). 


4.3.  Comparison  of  asymptotic  solutions  with  the 
numerical  result 


Fig.  12  presents  the  numerical  result,  together  with  asymp¬ 
totic  solutions  for  two  limits  of  x/k.  Eqs.  (60)  and  (63)  are  in 
good  agreement  with  the  numerical  result  for  nearly  all  ranges 
of  x.  In  most  regions,  except  at  the  start  of  channel,  the  large  x/k 
condition  can  be  applied,  which  produces  the  following  formula 
for  the  change  in  concentration  of  species  i: 


cl(x,  —1)  ~  —A 


.  x\  Wo/ h  +  0.5526 tan 0 

0.722  (  -  +  0.730 — — - 

kJ  Wo/  h  2  tan  0 


(64) 


5.  Conclusion 

A  two-dimensional  model  was  developed  to  imitate  channel 
flow  in  a  fuel  cell.  The  effect  of  the  inlet  conditions  and  the 
channel  geometry  on  the  concentration  distribution  was  inves¬ 
tigated  in  the  case  of  wall  suction  or  injection.  The  product  of 
the  Reynolds  number  and  the  Schmidt  number,  k,  was  useful  in 
manipulating  the  concentration  distribution  within  the  channel. 
For  k  >  300,  the  upper  part  of  the  channel  was  hardly  affected 
by  the  reaction  at  the  bottom.  The  draft  angle  of  the  channel 
was  shown  to  suppress  the  change  in  concentration  at  the  reac¬ 
tion  site.  The  performance  of  the  fuel  cell  for  different  values  of 
the  draft  angle  requires  further  study,  together  with  the  water- 
exhausting  capability  of  a  non-rectangular  channel.  For  different 
limits  of  x/k ,  two  asymptotic  solutions  were  found  that  showed 
good  agreement  with  the  numerical  result. 
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